Finite-temperature properties of Pb(Zri_j:Tia;)03 alloys from first principles 
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A first-principles-derived approach is developed to study finite-temperature properties of 
Pb(Zri_a;Tia;)03 (PZT) solid solutions near the morphotropic phase boundary (MPB). Structural 
and piezoelectric predictions are in excellent agreement with experimental data and direct first- 
principles results. A low-temperature monoclinic phase is confirmed to exist, and is demonstrated 
to act as a bridge between the well-known tetragonal and rhombohedral phases delimiting the MPB. 
A successful explanation for the large piezoelectricity found in PZT ceramics is also provided. 
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Ferroelectric perovskite A(B'B")03 alloys are of grow- 
ing importance for a variety of device applications 
and are also of great current fundamental interest since 
little is known about the effects responsible for their 
anomalous properties. A good example of an A(B',B")03 
solid solution that is of both fundamental and technolog- 
ical importance is the Pb(Zri_a;Tia;)03 system. Usually 
denoted as PZT, this mixed-cation alloy is currently in 
widespread use in piezoelectric transducers and actua- 
tors |l| . Its phase diagram exhibits a morphotropic phase 
boundary (MPB) separating a region with a tetragonal 
ground state {x > 0.52) from a region with rhombohedral 
symmetry {x < 0.45) 

High piezoelectric response is experimentally found in 
ceramics of PZT around the MPB. The origins of this 
large piezoelectric response are unclear. On the one hand, 
semi-empirical simulations predict that the large exper- 
imental value of the c?33 piezoelectric coefficient results 
mainly from the large value of ^33 that a single-crystal 
PZT would exhibit On the other hand, recent first- 
principles calculations |^,^ have found that the ^33 coef- 
ficient of a tetragonal single crystal of Pb(Zro.5Tio.5)03 
are estimated to be three times smaller than the experi- 
mental value obtained for ceramics at low temperature. 

Furthermore, recent synchrotron x-ray powder diffrac- 
tion studies have revealed the existence of an unexpected 
low-temperature monoclinic phase of PZT at a;=0.48 0, 
which implies that the phase diagram of PZT is more 
complex than previously thought. This monoclinic phase 
may act as a second-order transitional bridge between the 
tetragonal phase, for which the electrical polarization P 
lies along the pseudo-cubic [001] direction, and the rhom- 
bohedral phase, for which P is along the pseudo-cubic 
[111] direction. If this is indeed the case, the polariza- 
tion of the monoclinic phase continuously rotates as the 
composition x decreases in the MPB region Such a 
continuous rotation has yet to be observed. 



Obviously, accurate simulations are needed to under- 
stand the properties of perovskite alloys in general, and 
of PZT in particular. Since the beginning of the present 
decade, first-principles methods have emerged as a pow- 
erful tool for investigating properties of ferroelectric sys- 
tems theoretically (see and references therein). 
However, these methods are essentially restricted to the 
study of the zero-temperature properties of small cells, 
while accurate and interesting predictions of alloy prop- 
erties would require calculations on much larger cells at 
finite temperature. Ideally one desires a computational 
scheme with the capability of predicting the properties 
of "real" perovskite alloy systems at finite temperature, 
with the accuracy of the first-principles methods. 

The purpose of this letter is to demonstrate that it 
is possible to develop such a scheme, and to apply it 
to study the finite-temperature behavior of PZT in the 
vicinity of the MPB. Remarkably, we find that the exis- 
tence of an intermediate monoclinic phase emerges natu- 
rally from this approach. Moreover, the theory provides 
a novel and successful explanation for the large piezoelec- 
tric response of PZT near the MPB, thereby explaining 
and resolving the previous theoretical difficulties in ob- 
taining agreement with the known experimental values 
of the piezoelectric coefficients. 

Our scheme is based on the construction of an effective 
Hamiltonian from first-principles calculations. A ferro- 
electric effective Hamiltonian |l0|] must include the fer- 
roelectric local soft mode and the strain variables, since 
ferroelectric transitions are accompanied by a softening of 
a polar phonon mode and by the appearance of a strain. 
An alloy effective Hamiltonian must also include the com- 
positional degrees of freedom. We propose to incorporate 
all such degrees of freedom by writing the total energy E 
as a sum of two energies, 

£:({uj,{v,},77i/,{crj}) = £'avo({uJ,{v,},77ff) 
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+ £;ioc({uJ,{v,},{a,}) , (1) 

where is the local soft mode in unit cell i; {v^} are the 
dimensionless local displacements which are related to 
the inhomogeneous strain variables inside each cell [ p^ ; 
r]H is the homogeneous strain tensor; and the {(Jj} char- 
acterize the atomic configuration of the alloy. That is, 
aj—+l or —1 corresponds to the presence of a B' or B" 
atom, respectively, at lattice site j of the A(B'i_j;B"a;)03 
alloy. The energy E'avo depends only on the soft mode and 
strain variables. The {cj} parameters are thus incorpo- 
rated into the second energy term -Eioc- 

For i?avo, we generalize the analytical expression suc- 
cessfully used in Ref. |]lO| for simple ABO3 systems to 
the case of an alloy. This generalization simply con- 
sists in using the virtual crystal approximation (VGA) 
^t| , i.e., in replacing A(B'i_j;,B"^)03 by a uniform but 
composition-dependent "virtual" ABO3 system. E^vc 
thus consists of five parts: a local- mode self-energy, a 
long-range dipole-dipole interaction, a short range inter- 
action between soft modes, an elastic energy, and an in- 
teraction between the local modes and local strain . 

For i?ioc, we propose an expression that includes (i) 
the onsite effect of alloying on the self-energy up to the 
fourth order in the local-mode vector ; and (ii) intersite 
contributions which are linear in and in v^: 

i;ioc({u.},{v,}, {cTj}) = 

i 

+ f^j* • u» + fji • V,)] , (2) 

where the sums over i and j run over unit cells and mixed 
sublattice sites, respectively. Aa(tTi) and A'^iui) charac- 
terize the onsite contribution of alloying, while Qj,iio'j) 
and Rj^i{crj) are related to alloying-induced intersite in- 
teractions. Here eji is a unit vector joining the site j to 
the center of the soft-mode vector u^, and tji is a unit 
vector joining the site j to the origin of the displacement 
Vi. In principle, terms involving higher powers of and 
Vi might be included to improve the quality of the ex- 
pansion, but as will be shown below, we find this level 
of expansion to give a very good account of experimental 
findings. We also find that Qj,i{(Jj) and Rj^i{aj) decrease 
rapidly as the distance between i and j increases. As a re- 
sult, we included contributions up to the third neighbors 
for Qj,i{(yj), and over the first-neighbor shell for Rj_i{(jj). 

All the parameters of Eqs. (1) and (2) are derived from 
first principles. The 18 parameters of i?avc (see Table II 
of Ref. |lO) ) are determined by fitting the results of first- 
principles VGA calculations. On the other hand, Aa(cri), 
^7(ci); s-iid Rj^i{uj) are derived by performing 

first-principles calculations in which a true atom [e.g., Ti 
or Zr in Pb(Zr,Ti)03] is surrounded by VGA atoms. The 
first-principles method used in the present study is the 
plane-wave ultrasoft-pseudopotential method within 
the local-density approximation (LDA) [^. The VGA 
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FIG. 1. Largest, middle, and smallest average cartesian 
coordinates u\, U2 and U3 of the local-mode vector as 
a function of temperature in disordered single crystals of 
Pb(Zro. 5X10.5)03. Solid lines denote inclusion, while dashed 
li nes denote neglect, of Eioc in Eq. (1). For clarity, U2 and 
U3 are not shown for the latter (VGA) case; they are nearly 
zero at all temperatures. 

approach adopted averages the B' and B" pseudopoten- 
tials, and is the one of Ref. Q. 

Once our effective Hamiltonian is fully specified, the 
total energy of Eq. (1) is used in Monte-Gar lo simula- 
tions to compute finite-temperature properties of ferro- 
electric alloys. We typically use a 12 x 12 x 12 super- 
cell (8640 atoms), since this choice yields well-converged 
results. The {ctj} variables are chosen randomly in or- 
der to mimic maximal compositional disorder - consis- 
tent with experimental reality [[Tst- and are kept fixed 
during the Monte-Garlo simulations. We find that aver- 
aging our results over a couple of different realizations 
of the disorder leads to well-converged statistical prop- 
erties. The outputs of the Monte-Garlo procedure are 
the local mode vectors u (directly related to the electri- 
cal polarization) , and the homogeneous strain tensor rjn . 
We use the correlation-function approach of Ref. ||l^ to 
derive the piezoelectric response from these Monte-carlo 
simulations. Up to 10^ Monte-Garlo sweeps are first per- 
formed to equilibrate the system, and then 2x 10^ sweeps 
are used to get the various statistical averages. The tem- 
perature is decreased in small steps. 

Figure 1 shows the largest, middle and smallest carte- 
sian coordinates (wi, U2 and U3) of the supercell average 
of the local mode vectors in Pb(Zro.5Tio.5)03 as a func- 
tion of the temperature, as predicted by our approach 
described by Eqs. (1) and (2). Each coordinate is close 
to zero at high temperature, characterizing a paraelectric 
cubic phase. As the system is cooled down, ui drastically 
increases while U2 and U3 remain nearly null. This indi- 
cates a transition to a ferroelectric tetragonal phase, con- 
sistent with measurements [|| . We predict that the spon- 
taneous polarization reaches 0.79 G/m^ at very low tem- 
perature, which compares well with the first-principles 
resuhs of 0.70 and 0.74 G/ve? [||. The tetragonal ax- 
ial ratio c/a ranges from 1 close to the transition region 
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FIG. 2. Piezoelectric coefficients as a function of tempera- 
ture in disordered (4mm) Pb(Zro.5Tio.5)03. Solid lines denote 
inclusion, while dashed lines denote neglect, of -Eioc in Eq. (1). 
Temperature has been rescaled to fit the experimental value 
of the Curie temperature. Statistical errors are estimated to 
be ~10% of the values displayed. 



to 1.02 for lower temperature. This is in good agree- 
ment both with the experimental value of 1.02 - 1.025 
1^,0 obtained for disordered samples, and with the first- 
principles result of 1.03 obtained for an ordered alloy 
Figure 1 also shows the predictions of the VGA alloy 
theory, corresponding to the neglect of Eioc in Eq. (1). 
Interestingly, one sees that i^ioc has no effect on the phase 
transition sequence, which is consistent with recent find- 
ings that the VGA appro ach ca n reproduce some struc- 
tural properties of PZT 



17[ . Whether or not Eioc 



is included in the total energy, we find a Gurie temper- 
ature Tc that is higher than the experimental value of 
640 K |l^ . This difhculty of reproducing Tc is a general 
feature of the effective-Hamiltonian approach ||l^,^,^ , 
and may be due to higher perturbative terms neglected 
in the analytical expression for the total energy. In order 
to compare our results with experimental data, we will 
henceforth rescale our temperature as in Ref. so that 
the theoretical Tc is forced to match the experimental 
one. 

Figure 2 shows the piezoelectric coefficients predicted 
for a tetragonal single crystal of Pb(Zro.5Tio.5)03 as a 
function of the rescaled temperature, when neglecting or 
incorporating i^ioc in Eq. (1). The independent coeffi- 
cients for the 4mm point group are ^33, ^31 and di^. 
One can notice that inclusion of ii^ioc has only a small 
effect on c?3i and ^33 : ^31 is rather small for any temper- 
ature except near the transition, and ^33 is around 50 - 
55 pC/N at room temperature in both simulations. Us- 
ing the experimental values of the elastic compliances ||] 
to compute 633 from our calculated d^i and o?33, we find 
that our alloy effective hamiltonian leads to an 633 of 4.3 
G/m^, while neglecting i^ioc in Eq. (1) yields a similar 633 
of 3.8 G/m^ at low temperature. Both predictions agree 
well with the first-principles results ranging between 3.4 
and 4.8 C/m^ |||,|jl4|, confirming that the VGA can re- 
produce the 633 coefficient of PZT 
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FIG. 3. Similar to Fig. (1), but computed as a function of 
composition x. The temperature of the simulation is 50 K, 
corresponding to a rescaled experimental value of 30 K (see 
text). 



Figure 2 also demonstrates that incorporating E'loc in 
the total energy leads to a large enhancement of the dis 
coefficient, which is consistent with recent measurements 
revealing that the piezoelectric elongation of the tetrag- 
onal unit cell of PZT does not occur along the polar di- 
rection 1^ . This enhancement is highly relevant for the 
piezoelectric response (^33 in ceramic samples, denoted 
dss^c, which involves an average of the form 



^33, c 



/ [(^31 + di5) sin^ 6 + ^33 cos^ 0] sin 6 cos dO 
Jo 

(3) 

over the single-crystal coefficients. The true alloy ap- 
proach of Eqs. (l)-(2) leads to a d^s^c of 163 pG/N at 
room temperature, in excellent agreement with the ex- 
perimental value of 170 pG/N pp^]. On the other hand, 
neglecting E'loc leads to a smaller dss.c of 90 pG/N. This 
difference clearly demonstrates the necessity of incorpo- 
rating the local alloying effect into the total energy for 
understanding the large piezoelectric response of PZT ce- 
ramics near the MPB. 

We now use our alloy effective Hamiltonian to inves- 
tigate the low-temperature phases of Pb(Zri_a;Tia;)03 
near the MPB. We choose a constant temperature of 50 
K in the Monte-Garlo simulations, and vary the over- 
all composition of the solid solution. This compositional 
variation affects two quantities: (i) the populations of 
aj equal to -|-1 or —1, and (ii) the alloy-related pa- 
rameters. For the latter, only the parameters entering 
the local-mode self-energy of i?ave fl^l; ^-^id the Aa{ai) 
and Aj{ai) in Eq. (2), are allowed to be composition- 
dependent. This composition-dependence is assumed to 
be linear, and is determined by performing first-principles 
simulations on two different compositional cells. Such a 
linear composition-dependence approach is only realistic 
when exploring a narrow range of compositions, as done 
in the present study. 

Figure 3 shows that the local mode, and hence the po- 
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larization, is parallel to the pseudo-cubic [001] direction 
for Ti compositions larger than 50%, which is consistent 
with a tetragonal phase. For compositions lower than 
47%, the polarization becomes parallel to the pseudo- 
cubic [111] direction, indicating the "high temperature" 
rhombohedral phase [^|J2^]. The most interesting fea- 
ture of Fig. 3 is the behavior of the local mode for the 
compositional range between 47.5% and 49.5%: as x 
decreases, ui decreases, while U2 and M3 increase and 
remain nearly equal to each other. This behavior is 
characteristic of an intermediate phase that is neither 
tetragonal nor rhombohedral. The strain tensor given by 
our simulations indicates that this intermediate phase is 
the monoclinic phase experimentally found by Noheda 
et al. We further predict that the monoclinic phase for 
X ~ 48% can be characterized by an angle of 90.7° and 
by lattice vectors — ao(— 1.005,— 1.005,— 0.009), = 
ao(1.002, -1.002,0.000), and c„ = ao(0.004,0.004,1.018), 
where ao is a cubic lattice constant. All these predictions 
are in excellent quantitative agreement with the experi- 
mental results of Ref. [Q] . Figure 3 clearly demonstrates 
that the monoclinic phase acts as a bridge between the 
rhombohedral and tetragonal phases, as indicated by the 
continuous rotation of the polarization as a function of 
composition. Our computational scheme is also able to 
reproduce the compositional range narrowing of the mon- 
oclinic phase observed when increasing the temperature 
It should be noted that a VGA-only calculation 
(i.e., neglecting i?ioc) does not reveal a monoclinic phase. 
This finding demonstrates once again the need for incor- 
porating the local effect of alloying into the total energy 
to study subtle effects. 

In summary, we have developed a first-principles de- 
rived computational scheme to study finite-temperature 
properties of Pb(Zri_j;Ti2;)03 solid solutions near the 
MPB as a function of composition and temperature. We 
find that there is a low-temperature monoclinic phase 
acting as a bridge between the rhombohedral phase, ex- 
isting for X < 0.47, and the tetragonal phase, occurring 
for X larger than 0.50. The predicted structural data are 
in very good agreement with measurements, as well as 
with direct first-principles calculations. The use of this 
approach also provides an explanation for the large ex- 
perimental value of (i33 in tetragonal ceramics of PZT 
near the MPB |2^. This large piezoelectricity is simply 
due to the very large value of the c?i5 coefficient predicted 
to occur in the single crystal. 
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